Tutorial: Running populations with binary_c-python (Monte-carlo sampling)

WARNING: the Monte-Carlo sampling functionality is still in beta-stage and is not fully reliable. This notebook is not finished yet either.

This notebook will show you how to evolve a population of stars through a Monte-Carlo based sampling method.

In the population evolution notebook we explained how to set up a Population object and evolve all the systems with a grid-based sampling method. If you are not familiar yet with the grid-based approach, we recommend studying that notebook first. Here we will use a Monte-Carlo based sampling approach.

This approach samples systems based on the probability of occurence according to their probability distributions. We thus do not assign a probability to the system, but systems with properties that have a higher probability of occurrence will be sampled more often than others. This is the most standard Monte-Carlo based sampling approach. More advanced features will be introduced later, like Latin-Hypercube sampling and Adaptive Importance sampling.

import os

from binarycpython.utils.custom_logging_functions import temp_dir
from binarycpython import Population

TMP_DIR = temp_dir("notebooks", "notebook_monte_carlo", clean_path=True)

# help(Population) # Uncomment to see the public functions of this object

We control which quantities we sample again by using sampling_variables like in the population evolution notebook. There are some important things to take care of here, like calculating the multiplicity at the correct point, and introducing branching_points at the correct variables.

Let’s look at an example of a single-star population, where we only sample the quantity M_1.

We will first set up the population with the correct settings.

# Create population object
example_pop = Population(tmp_dir=TMP_DIR, verbosity=1)

# Configure population object
example_pop.set(
    # binary_c physics options
    multiplicity=1,
    evolution_type="monte_carlo",
    monte_carlo_mass_threshold=10,
    _actually_evolve_system=False, # For now, let's disable actually evolving the systems with binary_c
    num_cores=1, # Lets set the number of cores to 1 for now
)

adding: evolution_type=monte_carlo to grid_options
adding: monte_carlo_mass_threshold=10 to grid_options
adding: _actually_evolve_system=False to grid_options
adding: num_cores=1 to grid_options

resolution = {'M_1': 10}

# add sampling variable for M_1
example_pop.add_sampling_variable(
    name="M_1",
    longname="Primary mass",
    valuerange=[0.06, 300],
    samplerfunc="10**self.const_linear(np.log10(0.10001), np.log10(300), {})".format(resolution["M_1"]),
    probdist="self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 301, -1.3, -2.3, -2.3)",
    bottomcode="multiplicity = self._sample_multiplicity(system_dict); system_dict['multiplicity'] = multiplicity;",
    dphasevol="dlnm1",
    parameter_name="M_1",
    branchpoint=1,
    branchcode="multiplicity == 1"
)

example_pop.evolve()

Warning: No parse function set. Make sure you intended to do this.
Creating and loading custom logging functionality
[2023-01-25 23:51:02,029 DEBUG    MainProcess] --- setting up the system_queue_filler now
Write sampling functions code to /tmp/binary_c_python-david/notebooks/notebook_monte_carlo/binary_c_monte_carlo_sampling_sampling_functions_bd5e5d94033f45cb9e26c1cd12ab7e9e.py
Load monte-carlo sampling functions from /tmp/binary_c_python-david/notebooks/notebook_monte_carlo/binary_c_monte_carlo_sampling_sampling_functions_bd5e5d94033f45cb9e26c1cd12ab7e9e.py
Writing sampling variable M_1 to monte-carlo sampling generator
Write grid code to /tmp/binary_c_python-david/notebooks/notebook_monte_carlo/binary_c_monte_carlo_sampling_generator_bd5e5d94033f45cb9e26c1cd12ab7e9e.py
Load monte-carlo generator function from /tmp/binary_c_python-david/notebooks/notebook_monte_carlo/binary_c_monte_carlo_sampling_generator_bd5e5d94033f45cb9e26c1cd12ab7e9e.py
Grid code loaded
[2023-01-25 23:51:02,055 DEBUG    Process-28] --- Setting up processor: process-0
[2023-01-25 23:51:02,080 DEBUG    MainProcess] --- Signalling processes to stop
Do join of subprocesses ...[2023-01-25 23:51:02,082 DEBUG    Process-28] --- Process-0 is finishing.

process 0 free memory and return
****************************************************
*                Process 0 finished:               *
*  generator started at 2023-01-25T23:51:02.030944 *
* generator finished at 2023-01-25T23:51:02.114483 *
*                   total: 0.08s                   *
*           of which 0.00s with binary_c           *
*                  Ran 33 systems                  *
*          with a total probability of 33          *
*         This thread had 0 failing systems        *
*       with a total failed probability of 0       *
*   Skipped a total of 0 zero-probability systems  *
*                                                  *
****************************************************

process 0 queue put output_dict
[2023-01-25 23:51:02,124 DEBUG    Process-28] --- Process-0 is finished.
process 0 return
Joined all subprocesses.
**********************************************************
*  Population-bd5e5d94033f45cb9e26c1cd12ab7e9e finished! *
*              The total probability is 33.              *
*  It took a total of 0.61s to run 33 systems on 1 cores *
*                   = 0.61s of CPU time.                 *
*              Maximum memory use 178.703 MB             *
**********************************************************

No failed systems were found in this run.
Do analytics
Added analytics to metadata

{'population_id': 'bd5e5d94033f45cb9e26c1cd12ab7e9e',
 'evolution_type': 'monte_carlo',
 'failed_count': 0,
 'failed_prob': 0,
 'failed_systems_error_codes': [],
 'errors_exceeded': False,
 'errors_found': False,
 'total_probability': 33,
 'total_count': 33,
 'start_timestamp': 1674690661.9999769,
 'end_timestamp': 1674690662.6125712,
 'time_elapsed': 0.6125943660736084,
 'total_mass_run': 10.188142293973625,
 'total_probability_weighted_mass_run': 10.188142293973625,
 'zero_prob_stars_skipped': 0}